An efficient computation of geometric entanglement for two-dimensional quantum lattice systems 
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The geometric entanglement per lattice site, as a holistic measure of the multipartite entanglement, serves 
as a universal marker to detect quantum phase transitions in quantum many-body systems. However, it is very 
difficult to compute the geometric entanglement due to the fact that it involves a complicated optimization 
over all the possible separable states. In this paper, we propose a systematic method to efficiently compute 
the geometric entanglement per lattice site for quantum many-body lattice systems in two spatial dimensions 
in the context of a newly-developed tensor network algorithm based on an infinite projected entangled pair 
state representation. It is tested for quantum Ising model in a transverse magnetic field and anisotropic spin 
1 12 anti-ferromagnetic XYX model in an external magnetic field on an infinite-size square lattice. In addition, 
the geometric entanglement per lattice site is able to detect the so-called factorizing field. Our results are in a 
quantitative agreement with Quantum Monte Carlo simulations. 
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Introduction. The geometrical entanglement is a measure of 
the multipartite entanglement present in a quantum state wave 
function. It quantifies the distance between a given quantum 
state wave function and the closest separable state [lLl2D- Re- 
markably, it serves as an alternative marker J2-0] to locate 
critical points for quantum many-body lattice systems under- 
going quantum phase transitions (QPTs) In addition, 
an intriguing connection to both the Renormalization Group 
and Conformal Field Theory has been unveiled for quantum 
many-body lattice systems in one spatial dimension [7-9]. 
Furthermore, recent numerical simulations [8] established a 
universal finite-size correction to the geometric entanglement 
for the critical XXZ and transverse quantum Ising chains, 
which in turn is related with the celebrated Affleck-Ludwig 
g factors lilOiJ]- Therefore, the geometrical entanglement of- 
fers a powerful tool to investigate quantum criticality in quan- 
tum many-body lattice systems. However, almost all studies 
have been restricted to quantum many-body lattice systems 
in one spatial dimension, with an exception Ultl . in which 
the geometrical entanglement is exploited to study QPTs for 
quantum many-body lattice systems in two spatial dimen- 
sions. This is mainly due to the fact that it is very difficult 
to compute the geometrical entanglement, because it involves 
a complicated optimization over all the possible separable 
states. 

On the other hand, significant progress has been made to 
develop efficient numerical algorithms to simulate quantum 
many -body lattice systems in the context of the tensor network 
representations Ill2l42lll . The algorithms have been success- 
fully exploited to compute the ground-state fidelity per lattice 
site 1122142611 . a universal marker to detect QPTs, for quantum 
many-body lattice systems. Indeed, the ground-state fidelity 
per lattice site is closely related to the geometrical entangle- 
ment. Therefore, it is natural to expect that there should be 
an efficient way to compute the geometrical entanglement in 
the context of the tensor network algorithms. Actually, this 
has been achieved for quantum many-body lattice systems 
with the periodic boundary conditions in one spatial dimen- 



sion ifioll in the context of the matrix product state representa- 
tion. 

In the present work, we propose a systematic method to ef- 
ficiently compute the geometric entanglement per lattice site 
for quantum many-body lattice systems in two spatial dimen- 
sions in the context of a newly-developed tensor network al- 
gorithm based on an infinite projected entangled pair state 
(iPEPS) representation. It is exploited to evaluate the geo- 
metric entanglement per lattice site for quantum Ising model 
in a transverse magnetic field and anisotropic spin 1/2 anti- 
ferromagnetic XYX model in an external magnetic field on an 
infinite-size square lattice, which enables us to identify their 
critical points. In addition, the geometric entanglement per 
lattice site is able to detect the so-called factorizing field. Our 
results are in a quantitative agreement with Quantum Monte 
Carlo simulations. 

The geometric entanglement per lattice site. For a pure 
quantum state \ifr) with N parties, the geometric entanglement, 
as a global measure of the multipartite entanglement, quanti- 
fies the deviation from the closest separable state \<p). Math- 
ematically, the geometric entanglement [H |9[] E(\i//)) for an 

-partite quantum state \i/r) is expressed as: 



£(|-A» = -io g2 A;L 
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where A max is the maximum fidelity between and all the 
possible separable and normalized state \<j>), with 
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Then, the geometric entanglement per party £#(|0)) is defined 
as: 

S N (\ifr)) = N-'Em). (3) 
It corresponds to the maximum fidelity per party /i™ x , where 

= V^W, (4a) 

or, equivalently, 

£ w (| ( A»--log 2 (^ ax ) 2 . (4b) 



2 



For our purpose, we shall consider a quantum many-body 
system on an infinite-size lattice in two spatial dimensions, 
which undergoes a QPT at a critical point in the thermody- 
namic limit. In this situation, each lattice site constitutes a 
party, thus the geometric entanglement per party becomes the 
geometric entanglement per lattice site, which is well defined 
even in the thermodynamic limit (N — > oo), since the contri- 
bution to fidelity from each party (site) is multiplicative. 

The infinite projected entangled pair state algorithm. Our 
aim is to compute the geometric entanglement per lattice site 
for a quantum many-body lattice system on an infinite-size 
square lattice in the context of the iPEPS algorithm [ 14]. Sup- 
pose we consider a system characterized by a translation- 
invariant Hamiltonian H with the nearest-neighbor interac- 
tions: H - Yi<ij> h<ij>, with /?</)> being the Hamiltonian den- 
sity. Assume that a quantum wave function \iff) is translation- 
invariant under two-site shifts, then one only needs two five- 
index tensors Af, and Bf, to express the iPEPS representa- 
tion. Here, each tensor is labeled by one physical index s and 
four bond indices I, r, u and d, as shown in FigJTJi). Note that 
the physical index s runs over 1, • • • , d, and each bond index 
takes 1 , • • • , D, with d being the physical dimension, and D 
the bond dimension. Therefore, it is convenient to choose a 
2x2 plaquette as the unit cell (cf. FigJTJii)). The ground-state 
wave function is well approximated by \i/f T ), which is obtained 
by performing an imaginary time evolution [14] from an ini- 
tial state \ifr ), with \ifr T ) = e - HT |<A()>/lk^ T |(Ao>ll H, as long 
as r is large enough. 

A key ingredient of the iPEPS algorithm is to take advan- 
tage of the Trotter-Suzuki decomposition that allows to reduce 
the (imaginary) time evolution operator e~ HSr over a time slice 
St into the product of a series of two-site operators, where the 
imaginary time interval r is divided into M slices: r = M8t. 
Therefore, the original global optimization problem becomes 
a local two-site optimization problem. With an efficient con- 
traction scheme available to compute the effective environ- 

1411 . one is able 
Performing this proce- 
dure until the energy per lattice site converges, the ground- 
state wave function is produced in the iPEPS representation. 

Efficient computation of the geometric entanglement in the 
iPEPS representation. Once the iPEPS representation for the 
ground-state wave function is generated, we are ready to eval- 
uate the geometric entanglement per lattice site. First, we need 
to compute the fidelity between the ground-state wave func- 
tion and a separable state. The latter is represented in terms 
of one-index tensors A s and B s . To this end, we form a re- 
duced four-index tensor from the five-index tensor Af, 
and a one-index tensor A s , as depicted in Fig. [TJ(iii). As such, 
the fidelity is represented as a tensor network in terms of the 
reduced tensors a\ ru d and bi n ,d (cf. Fig. Q](iv)). The tensor 
network may be contracted as follows. First, form the one- 
dimensional transfer matrix E\, consisting of two consecu- 
tive rows of the tensors in the checkerboard tensor network. 
This is highlighted in Fig.[T](vi) with two dash lines. Second, 
compute the dominant eigenvectors of the transfer matrix E\, 
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FIG. 1: (color online) (i) A five-index tensor A? , labeled by one 
physical index s and four bond indices /, r, u and d. (ii) The iPEPS 
representation of a wave function in two-dimensional square lattice. 
Copies of the tensors A s lmd and B*, are connected through four types 
of bonds, (iii) A one-index tensor A 1 labeled by one physical index s. 
(iv) The iPEPS representation of a separable state in two-dimensional 
square lattice, (v) A reduced four-index tensor «/,,«/ from a five-index 
tensor Af . and a one-index A", (vi) The tensor network represen- 
tation for the fidelity between a quantum wave function (describe by 
Af , and B s lmd ) and a separable state (described by A 5 and B s ), con- 
sisting of the reduced tensors a\ rui and b lrud . 



corresponding to the dominant eigenvalue. This can be done, 
following a procedure described in Ref. [lg]. Here, the dom- 
inant eigenvectors are represented in the infinite matrix prod- 
uct states. Third, choose the zero-dimensional transfer ma- 
trix Eq, and compute its dominant left and right eigenvectors, 
Vl and Vr. This may be achieved by means of the Lanczos 
method. In addition, one also needs to compute the norms 
of the ground-state wave function \ifr) and a separable state |0) 
from their iPEPS representations. Putting everything together, 
we are able to get the fidelity per unit cell between the ground 
state \iff) and a separable state \<p): 



A = 



(5) 



where tj^), t]^],/,) and t]^) are, respectively, the dominant 
eigenvalue of the zero-dimensional transfer matrix Eq for the 
iPEPS representation of <0|i/'), and (<p\(f>). Then we pro- 
ceed to compute the geometric entanglement per lattice site, 
which involves the optimization over all the separable states. 
For our purpose, we define F — A 2 . The optimization amounts 
to computing the logarithmic derivative of F with respect to 
A*, which is expressed as 



G = 



1 



d\nF 1 d-q m) _ 
dA* 77<0|<fr> dA* 77(0|0> dA* 



(6) 



Therefore, the problem reduces to the computation of G in 
the context of the tensor network representation. First, note 
that a pictorial representation of the derivative dai ru d/dA s * of 
the four-index tensor ai ru d with respect to A s * is shown in 
Fig. |2] (ii), which is nothing but the five-index tensor Af 
Similarly, we may define the derivative of the four-index ten- 
sor birud with respect to B s *. Then, we are able to represent 
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FIG. 2: (color online) The gradient of the fidelity between a given 
ground-state wave function \tf/) and a separable state \(p) in the iPEPS 
representation, (i): The zero-dimensional transfer matrix Eq and its 
dominant eigenvectors V L and V R . Here, the infinite matrix prod- 
uct state representation of the dominant eigenvectors for the one- 
dimensional transfer matrix E\ follows from Ref. 1 16], and V L and V R 
may be evaluated from the Lanczos method. The contraction of the 
entire tensor network is the dominant eigenvalue rj^) of the zero- 
dimensional transfer matrix Eq for (if>\if/). (ii): A half-filled square 
denotes a_, the derivative of the four-index tensor a h . ui with respect 
to i", which is nothing but the five-index tensor A? .. Similarly, 
we may define the derivative of the four-index tensor b lmd with 
respect to B s * . (iii) and (iv): The pictorial representation of the con- 
tributions to the derivative of rj^y with respect to A", with different 
relative positions between filled circles and half-filled squares. 



the contributions to the derivative of rj^) with respect to A 1 * 
in Fig. |2] (iii) and (iv). In our scheme, we update the real and 
imaginary parts of A s separately: 



and 



%{A S ) = %(A S ) + 6%(G) S 



3(A S ) = 3(A S ) + S3(G) S 



(7a) 



(7b) 



Here, 5 e [0, 1) is the step size in the parameter space, which 
is tuned to be decreasing during the optimization process. In 
addition, we have normalized the real and imaginary parts of 
the gradient G so that their respective largest entries are unity. 
The procedure to update the tensor B s is the same. If the fi- 
delity per unit cell converges, then the closest separable state 
|0) is achieved, thus the geometric entanglement per lattice 
site for the ground-state wave function follows. 

The models. As a test, we simulate two typical quantum 
many-body lattice systems on an infinite-size square lattice 
in two spatial dimensions. The first model is quantum Ising 
model in a transverse magnetic field. The Hamiltonian takes 
the form, 
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FIG. 3: (color online) The geometric entanglement per lattice site 
&oe(h) as a function of the transverse magnetic field strength h for 
quantum transverse Ising model on an infinite-size square lattice in 
two spatial dimensions. The data are presented for both the trun- 
cation dimension D = 2 and D = 3. The cusp in the S m (h) curve 
reflects a drastic change of the multipartite entanglement around a 
critical point: h c = 3.1 for D = 2 and h c = 3.08 for D = 3. Quantum 
Monte Carlo simulation indicates a critical point at /iqmc ~ 3.044. 



where S„ (a = x, z) are the spin- 1/2 Pauli operators at lattice 
site j, his a transverse magnetic field, and < i, j > runs over all 
the possible pairs of the nearest neighbors on a square lattice. 
Quantum Monte Carlo simulation predicts a critical point at 
h QMC ~ 3.044 S. 

The second model is quantum anti-ferromagnetic XYX 
model in an external magnetic field. The Hamiltonian is writ- 
ten as, 



H = J £ (sfS f + hySfsf + S®SW) + h ^ S M, (9) 



<>,]> 



where Sa(a = x,y, z) are the spin- 1/2 Pauli operators at lat- 
tice site j, J > is the exchange coupling, < i, j > runs over 
all the possible pairs of the nearest neighbors on a square lat- 
tice, and h is an external magnetic field. Note that A v < 1 
and A v > 1 correspond to the easy-plane and easy-axis behav- 
iors, respectively. Here, we focus on A v = 0.25. In this case, 
the model is known to undergo a continuous QPT in the same 
universality class as the transverse Ising model M28I1 . 

Simulation results. We plot the geometric entanglement per 
lattice site for quantum Ising model in a transverse field, with 
the field strength h as the control parameter, in Fig. [3] There 
is a cusp at h c = 3.10 for the truncation dimension D = 2 and 
h c = 3.08 for the truncation dimension D = 3, respectively, on 
the curve of the geometric entanglement as a function of the 
transverse field h. On both sides of the cusp, it is continuous. 
This implies that a continuous QPT occurs at h c . 

In Fig. |H the geometric entanglement per lattice site for 
quantum XYX model in an external magnetic field h, with 
A v = 0.25, is plotted. Here, the external magnetic field 
strength h is chosen as the control parameter. A cusp occurs as 
h varies across a critical point h c : h c = 3.49 for the truncation 
dimension D — 2 and h c = 3.485 for the truncation dimension 
D — 3, respectively. On both sides of the cusp, the geometric 
entanglement per lattice site is continuous, which implies that 



(8) 
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FIG. 4: (color online) The geometric entanglement per lattice site 
&co(h) as a function of the applied external magnetic field h for quan- 
tum XYX model on an infinite-size square lattice in two spatial di- 
mensions, with A v = 0.25. The data are presented for both D = 2 
and D = 3. The factorizing field, at which the geometric entangle- 
ment per lattice site vanishes, is indicated by an arrow labeled by 
hf around 3.162. Above the factorizing field hf, a sharp increase 
of £„(h) reflects a rapid increase of the multipartite entanglement 
around a critical point: h c = 3.49 for D = 2 and h c = 3.485 for 
D = 3. 



the model undergoes a continuous QPT at h c . 

Notice that, only a small shift is observed for h c , with 
increasing of the truncation dimension, for both models. 
This indicates that the iPEPS algorithm captures many-body 
physics for both models, with a small truncation dimension, as 
already observed before 11141 l23l 12611 . Our simulation results 
are in a quantitative agreement with Quantum Monte Carlo 
simulation for quantum Ising model in a transverse field 12711 
and for quantum XYX model in an external magnetic field h, 
with A v = 0.25 12811 . In addition, our simulation indicates that 
a factorizing field hf occurs at hf = and hf = 3.162, re- 
spectively, for quantum Ising model in a transverse field and 
for quantum XYX model in an external magnetic field h, with 
A,, = 0.25, thus reproducing the exact results. We stress that a 
factorizing field hf is identified fromfi M (/i/) = 0. 

Conclusions. In this paper, we have demonstrated how 
to efficiently compute the geometric entanglement per lat- 
tice site, by optimizing over all the possible separable states, 
in the context of the tensor network algorithm based on the 
iPEPS representation. The geometric entanglement per lattice 
site, as a holistic measure of the multipartite entanglement, 
serves as a universal marker to locate critical points underly- 
ing quantum many-body lattice systems. Our method is tested 
for both quantum Ising model in a transverse magnetic field 
and an anisotropic spin 1 /2 anti-ferromagnetic XYX model 
in an external magnetic field on an infinite-size square lattice, 
succeeded in identifying both the critical points and factor- 
izing fields. We expect that, with the developments of pow- 
erful tensor network algorithms, the geometric entanglement 
per lattice site adds a new route to explore quantum criticality 
for quantum many-body lattice systems in condensed matter 
physics. 

Acknowledgements. The work is partially supported by 



[1] H. Barnum and N. Linden, J. Phys. A: Math. Gen. 34, 6787 
(2001). 

[2] T.-C. Wei and P. M. Goldbart, Phys. Rev. A 68, 042307 (2003); 

T.-C. Wei, D. Das, S. Mukhopadyay, S. Vishveshwara, and P.M. 

Goldbart, Phys. Rev. A 71, 060305 (2005). 
[3] R. Orus and T.-C. Wei, Phys. Rev. B 82, 155120 (2010). 
[4] R. Orus, T.-C. Wei, and H.-H. Tu, arXiv: 1010.5029; R. Orus 

and T.-C. Wei, arXiv: 1006.5584. 
[5] S. Sachdev, Quantum Phase Transitions, Cambridge University 

Press, 1999, Cambridge. 
[6] X.-G. Wen, Quantum Field Theory of Many-Body Systems, Ox- 
ford University Press, 2004, Oxford. 
[7] A. Botero and B. Reznik, arXiv:0708.3391; R. Orus, Phys. Rev. 

Lett. 100, 130502 (2008). 
[8] Q.-Q. Shi, R. Orus, J. O. Fjasrestad, and H.-Q Zhou, New J. 

Phys. 12, 025008 (2010). 
[9] J.-M. Stephan, G. Misguich, and F. Alet, Phys. Rev. B 82, 

180406R (2010). 
[10] B.-Q. Hu, X.-J. Liu, J.-H. Liu, and H.-Q. Zhou, preprint. 
[11] C.-Y. Huang and F.-L. Lin, Phys. Rev. A 81, 032304 (2010). 
[12] F. Verstraete, D. Porras, and J.I. Cirac, Phys. Rev. Lett. 

93, 227205 (2004). F. Verstraete and J.I. Cirac, arXivxond- 

mat/0407066; V. Murg, F. Verstaete, and J.I. Cirac, Phys. Rev. 

A 75, 033605 (2007). 
[13] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003); G. Vidal, Phys. 

Rev. Lett. 93, 040502 (2004). G. Vidal, Phys. Rev. Lett. 98, 

070201 (2007). 

[14] J. Jordan, R. Orus, G. Vidal, F. Verstraete, and J. I. Cirac, Phys. 

Rev. Lett. 101, 250602 (2008). 
[15] B. Bauer, G. Vidal, and M. Troyer, J. Stat. Mech. (2009) 

P09006. 

[16] R. Orus and G. Vidal, Phys. Rev, B 78, 155117 (2009). 

[17] P. Pippan, S. R. White, and H. G. Evertz, Phys. Rev. B 81, 

081 103(R) (2010). 
[18] Q.-Q. Shi and H.-Q. Zhou, J. Phys. A: Math. Theor. 42 272002 

(2009). 

[19] J. Haegeman, J. I. Cirac, and T. J. Osborne, I. Pizorn, H. Ver- 

schelde, and F. Verstraete, arXiv: 1103.0936. 
[20] H. C. Jiang, Z. Y. Weng, and T. Xiang, Phys. Rev. Lett. 101, 

090603 (2008). 

[21] L. Wang, Y.-J. Kao, and A.W. Sandvik, Phys. Rev. E 83, 056703 
(2011). 

[22] H.-Q. Zhou and J.P. Barjaktarevic, J. Phys. A: Math. Theor. 41 
412001 (2008); H.-Q. Zhou, J.-H. Zhao, and B. Li, J. Phys. A: 
Math. Theor. 41 492002 (2008). 

[23] H.-Q. Zhou, R. Orus, and G. Vidal, Phys. Rev. Lett. 100, 
080602 (2008). 

[24] J.-H. Zhao, H.-L. Wang, B. Li, and H.-Q. Zhou, 
arXiv:0902.1669. Phys. Rev. E. 82, 061127 (2010). 

[25] H.-L. Wang, J.-H. Zhao, B. Li, and H.-Q. Zhou, 
arXiv:0902.1670. 

[26] B. Li, S.-H. Li, and H.-Q. Zhou, Phys. Rev. E 79, 060101(R) 
(2009). 



5 



[27] H.W.J. Blote and Y. Deng, Phys. Rev. E 66, 0661 10 (2002). 
[28] T. Roscilde et. al, Phys. Rev. Lett. 94, 147208 (2005). 



